/*
    NOTE: PSID_incscores_bothparents is created in the 
          build/PSID/code/Actual_vs_predictedincome
          folder. Run 1_PSID and 2_PSID do files first.
*/

clear 
set more off


   use "$PSID/output/PSID_incscores_bothparents.dta", clear

    levelsof birthyear_bin5, local(levels) 
    di "`levels'"

*--------------------------------------------------------*
* WORKING DAD SAMPLE
*--------------------------------------------------------*

    global workdad_incvars "log_mean_totfaminc_3years log_mean_fathertotfaminc_3years log_father_interpolated"

    preserve

      keep if analysis_samp_3yrs_workingdads==1 

    *--------------------------*
    * Create ranked measures
    *--------------------------*
        foreach var of global workdad_incvars  { 
        
            if "`var'" == "log_mean_totfaminc_3years" local rname0 "c_rankactual_wkdads_3yrs"
            if "`var'" == "log_mean_totfaminc_3years" local ysel0 "ysel_c_rankactual_wkdads_3yrs"
            
            if "`var'" == "log_mean_fathertotfaminc_3years" local rname0 "p_rankactual_wkdads_3yrs"
            if "`var'" == "log_mean_fathertotfaminc_3years" local ysel0 "ysel_p_rankactual_wkdads_3yrs"

            if "`var'" == "log_father_interpolated" local rname0 "p_rankpred_wkdads"
            if "`var'" == "log_father_interpolated" local ysel0 "ysel_p_rankpred_wkdads"
            
                egen `rname0' = xtile(`var') if inrange(birthyear,1946,1981), by(birthyear_bin5) nq(100) weight(weight1997)
                
                qui:gen `ysel0'=.
                replace `ysel0'= birthyear if `rname0'!=.

            }

        assert c_rankactual_wkdads_3yrs==. if log_mean_totfaminc_3years==. 
        assert p_rankactual_wkdads_3yrs==. if log_mean_fathertotfaminc_3years==.
        assert p_rankpred_wkdads==. if log_father_interpolated==.
        
        label var c_rankactual_wkdads_3yrs "Rank adult child, 3 yrs actual family income, working dad sample"
        label var p_rankactual_wkdads_3yrs "Rank dad, 3 yrs actual family income, working dad sample"
        label var p_rankpred_wkdads "Rank dad, predicted income, working dad sample"

    *--------------------*
    * REGRESSIONS
    *--------------------*
    /* First loop: actual dad ranked income 
       Second loop: predicted dad ranked income */

    foreach pvar in p_rankactual_wkdads_3yrs p_rankpred_wkdads {   
        if "`pvar'"=="p_rankactual_wkdads_3yrs" local name "rr_wkdad_actualinc"
        if "`pvar'"=="p_rankpred_wkdads" local name "rr_wkdad_predinc"

        matrix `name' = J(7,3,.)
        local i =0

        foreach y in `levels' {
            reg c_rankactual_wkdads_3yrs `pvar' if birthyear_bin5==`y' [w=weight1997], cluster(parent_id)

            local i=`i'+1
            matrix `name'[`i',1] = `y'
            matrix `name'[`i',2] = _b[`pvar']
            matrix `name'[`i',3] = _se[`pvar']
        }
    }

   restore

*---------------------------------------------------------------------------*
*---------------------------------------------------------------------------*

*------------------------------------*
* FIGURE
*------------------------------------*
         
    //Working dads, actual vs predicted (ranked) dad income
    clear 
    svmat rr_wkdad_actualinc
    ren rr_wkdad_actualinc1 bin
    gen l95_rr_wkdad_actualinc = rr_wkdad_actualinc2-1.96*rr_wkdad_actualinc3
    gen u95_rr_wkdad_actualinc = rr_wkdad_actualinc2+1.96*rr_wkdad_actualinc3
    sort bin
    tempfile rr_wkdad_actualinc
    save `rr_wkdad_actualinc'

    clear 
    svmat rr_wkdad_predinc
    ren rr_wkdad_predinc1 bin
    gen l95_rr_wkdad_predinc = rr_wkdad_predinc2-1.96*rr_wkdad_predinc3
    gen u95_rr_wkdad_predinc = rr_wkdad_predinc2+1.96*rr_wkdad_predinc3
    sort bin

    merge 1:1 bin using `rr_wkdad_actualinc'
    drop _merge

    gen bin2 = bin + 0.05

    preserve
        keep if bin<=3 //birth cohorts up to 1960
        #delimit ;
        twoway (connect rr_wkdad_actualinc2 bin, msymbol(circle_hollow) mcolor(pink*1.5) msize(medium) lcolor(pink*1.5) lwidth(0.25) lpat(solid)) (rcap u95_rr_wkdad_actualinc l95_rr_wkdad_actualinc bin, lcolor(pink*1.5) lwidth(0.25))
               (connect rr_wkdad_predinc2 bin2, msymbol(square) mcolor(lavender) msize(medium) lcolor(lavender) lwidth(0.25) lpat(solid)) (rcap u95_rr_wkdad_predinc l95_rr_wkdad_predinc bin2, lcolor(lavender) lwidth(0.25)),        
        graphregion(color(white)) legend(on size(small) ring(0) row(2) pos(5) order(1 "Actual father income" 3 "Predicted father income")) 
        xtitle(" " "Birth Year (binned)") ytitle("Rank-rank slope" " ", axis(1)) ylabel(0(0.1)0.6,nogrid axis(1)) xscale(range(0.5 3.5))
        xlabel(1 "1946-1950" 2 "1951-1955" 3 "1956-1960", labsize(small)  angle(45))  
        ;
        #delimit cr
        graph export "$Mydirectory2/appendix_d/rr_workingdads_PSID_upto1960.pdf", as(pdf) replace
    restore
